Liquid state properties from first principles DFT calculations: Static properties 
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In order to test the vibration-transit (V-T) theory of liquid dynamics, ab initio density functional 
theory (DFT) calculations of thermodynamic properties of Na and Cu are performed and compared 
with experimental data. The calculations are done for the crystal at T = and T m , and for the 
liquid at T m . The key theoretical quantities for crystal and liquid are the structural potential and 
the dynamical matrix, both as functions of volume. The theoretical equations are presented, as well 
as details of the DFT computations. The properties compared with experiment are the equilibrium 
volume, the isothermal bulk modulus, the internal energy and the entropy. The agreement of theory 
with experiment is uniformly good. Our primary conclusion is that the application of DFT to V-T 
theory is feasible, and the resulting liquid calculations achieve the same level of accuracy as does 
ab initio lattice dynamics for crystals. Moreover, given the well established reliability of DFT, the 
present results provide a significant confirmation of V-T theory itself. 

PACS numbers: 05.70.Ce, 64.10.+h 



I. INTRODUCTION 

The goal of our work is to investigate, and to improve 
where possible, the theoretical procedures for calculat- 
ing statistical mechanical properties of condensed matter 
systems. Here we shall focus on elemental crystals and 
liquids. For many such systems, ab initio density func- 
tional theory (DFT) provides highly accurate results for 
the groundstate energy as function of the nuclear posi- 
tions. This energy is the groundstate adiabatic potential, 
which appears in the nuclear motion Hamiltonian. For 
crystals, the nuclear motion Hamiltonian is prescribed by 
lattice dynamics theory [l|. Thermodynamic properties 
of elemental crystals, as calculated from DFT together 
with phonon statistical mechanics, can be nearly as ac- 
curate as the experimental measurements. For liquids, in 
the absence of a tractable nuclear motion Hamiltonian, 
statistical mechanical properties are calculated from ab 
initio molecular dynamics (MD). These calculations are 
based on DFT, and can be evaluated in the groundstate 
adiabatic approximation [2|, |3|. Again the results com- 
pare very well with experiment. Moreover, the studies 
reveal detailed characteristics of the electronic structure, 
for example for liquid Al [4J , for liquid Fe [9, Q , and for 
group III B - VI B elemental liquids J?]]- 

In recent years, vibration-transit (V-T) theory has 
been under development to provide a tractable approx- 
imate Hamiltonian for monatomic liquids @, Q • In this 
theory the nuclear motion is composed of two parts, the 
many-body vibrational motion in one potential energy 
valley, plus transits, which carry the system from val- 
ley to valley. Transits proceed at a high rate through- 
out the liquid, and are responsible for diffusion. The 
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vibrational motion is tractable and is subject to ab ini- 
tio evaluation. Closed-form equations are available for 
the dominant structural and vibrational contributions to 
thermodynamic functions. The transit contribution is 
complicated but small, and is treated by a model. The 
question we ask is simple and direct: If we apply DFT 
in the adiabatic approximation to V-T theory, how do 
the calculated thermodynamic properties compare with 
experiment? Our purpose here is to provide an initial 
answer to this question. 

The study is done for Na and Cu. The properties we 
calculate are the equilibrium volume, the isothermal bulk 
modulus, and the internal energy and entropy. Compar- 
ison of theory and experiment is done for the crystal at 
T = and at the melting temperature T m , and for the 
liquid at T m . The crystal tests are made to establish a 
fiducial for the theoretical accuracy. Both Na and Cu 
have a nearly-free-electron density of states in the vicin- 
ity of the Fermi energy. Hence the electronic excitation 
contribution to thermodynamic properties is very small 
and may be calculated from free electron theory. In this 
way, the ultimate comparison of theory and experiment 
is not significantly contaminated by error from electronic 
excitation. On the other hand, while Na has a rigid core, 
with only one valence electron that deforms as the nu- 
cleus moves, both the valence s-electron and the filled 
d-shell deform as the nucleus moves in Cu. This strong 
difference in the groundstate adiabatic potential adds di- 
mension to the present study. While this study is limited 
to two liquid metals, extensive analysis of experimental 
data reveal a common behavior of most elemental liquids 
[8J, and Na and Cu are representative of this common 
behavior. 

Before completing the work reported here, two pre- 
liminary results were required. First, an efficient DFT 
quench procedure for locating the many-body potential 
energy minima was developed [lOL |ll|. Second, an accu- 



rate predictive model for the transit contribution to ther- 
modynamic functions was constructed 12]. The culmina- 
tion of the complete project is reported here. It must be 
observed that the present work is not intended to replace 
ab initio MD calculations; indeed the two methods are 
quite complementary, as discussed in the final section. 

In Sec. |nl we outline the theory and explain various 
issues relevant to real- world calculations. In Sec. IIII1 the 
DFT calculations are described, including quench proce- 
dures and the total energy for crystal and liquid struc- 
tures. (Miscellaneous details related to these sections 
are collected in the appendices.) Results are presented 
and discussed in Sec. [TV] Intermediate theoretical results 
show the relative importance of the separate theoretical 
contributions to internal energy and entropy. Theory and 
experiment are tabulated and compared at the precision 
of the experimental accuracy. Conclusions are summa- 
rized in Sec.[VJ The primary conclusion regards the over- 
all accuracy of the present calculations of thermodynamic 
properties of liquid Na and Cu. A secondary conclusion 
summarizes the verification of V-T theory which is pro- 
vided by the present calculations. 



II. THEORETICAL FORMULATION 

The condensed matter potential energy surface is com- 
posed of intersecting many-atom potential energy valleys. 
Each valley makes a contribution to the partition func- 
tion. To approximate the single valley partition function, 
the valley potential is harmonically extended to infin- 
ity. The partition function is then simple, at the cost of 
neglecting anharmonicity and valley- valley intersections. 
We summarize the harmonic single valley statistical me- 
chanical formulas in Appendix[AJ In this section we show 
why these formulas are needed for the present work. 

For liquids, the starting point of V-T theory is a hy- 
pothesis about the nature of the many-body potential 
energy valleys which underlie the nuclear motion. These 
valleys are divided into two classes, random and sym- 
metric. In the thermodynamic limit, the random valleys 
are supposed to dominate the liquid statistical mechan- 
ics, and are also supposed to be macroscopically uniform. 
Macroscopic uniformity means that the statistical me- 
chanical average of any macroscopic dynamical variable 
is the same for all random valleys. Therefore the vibra- 
tional contribution to a thermodynamic function can be 
calculated from a single random valley harmonically ex- 
tended to infinity. The single-valley vibrational motion 
is supposed to be interspersed with transits, which carry 
the system from valley to valley. Hence there are two sep- 
arate components of the nuclear motion, vibrations and 
transits. With a superscript I to represent the liquid, the 
total free energy F l (V,T) is written 

F l (V, T) = & (V) + Fl ih (V, T) + FUV, T) + F^(V, T). 

(1) 
Here &q(V) is the system potential energy at the ran- 
dom valley structure, and F l vih {V,T) and Fl r (V,T) ex- 



press respectively the nuclear motion contribution from 
vibrations and transits. The final term, F^(V,T), repre- 
sents electronic excitations; it is added here because Na 
and Cu, the materials we study, are metals. However, 
F^(V,T) is very small, and free-electron theory in the 
leading Sommerfeld expansion provides sufficient accu- 
racy. 

The corresponding internal energy U (V, T) and en- 
tropy S l (V,T) are 



U l (V,T) -- 


= & (V) + UU(V,T) 






+UUV,T) + UUV,T), 


(2) 


S l {V,T) -- 


= S l vih (V,T) + S l u .(V,T) + S l cl (V,T). 


(3) 



At T > T m , the melting temperature, the high-T expan- 
sions Eqs. (|A10[) and (JA11I) are valid, hence the primary 
potential energy parameters are §b(V) for the internal 
energy and l Q (V) for the entropy. l Q (V) is the char- 
acteristic temperature related to the log moment of the 
vibrational spectrum (see Appendix [A)) . 

A new challenge, unique to the liquid, is to find the 
structures corresponding to random valleys so the char- 
acteristic functions may be calculated. The technique 
we've developed to do so [ljj,[ll[ exploits their numerical 
dominance; we start with computer generated stochas- 
tic configurations, in which the nuclei are distributed 
uniformly over the system volume, within a constraint 
limiting the closeness of approach of any pair. As we 
demonstrate in [Tfl Hl| . quenching from such a config- 
uration lands the system in a random valley with high 
likelihood. 

Once the structure is found, the system potential is 
corrected to the thermodynamic zero to produce ^ l a (V) 
(see Appendix |B|) . The dynamical matrix is the mass- 
weighted curvature tensor evaluated at the structure. 
This is calculated by a finite-difference approach in which 
each individual nucleus is displaced in all three Cartesian 
directions and the forces on all nuclei are computed. The 
eigenvalues are Muj\ for A = 1, . . . , 37V, where M is the 
nuclear mass and ui\ are the normal mode vibrational 
frequencies. Here, to calculate thermodynamic functions, 
only the eigenvalues are needed. However, the eigenvec- 
tors are also important, as they are needed to calculate 
fluctuations and time correlation functions [13| . 

Although we will perform these calculations in Sec. IIIII 
using periodic boundary conditions, that does not imply 
that the liquid vibrations correspond to those of a crystal 
with a large unit cell. The liquid vibrational modes are 
fundamentally different from those of a crystal. Since a 
crystal is periodic in space, periodic boundary conditions 
on a crystal unit cell or supercell merely express the in- 
finite extension of the crystal. This is the infinite lattice 
model [lj, an d entails no error. A liquid is represented 
by a random structure, which has no spatial periodicity, 
and calculations for a finite system contain surface errors. 
Such errors are minimized by the application of periodic 
boundary conditions, but the resulting spatial periodicity 
is not physically correct or meaningful for the liquid. 



TABLE I: Transit contributions to energy and entropy at 



melt. [Tj] 








Element 


T m /6ti 


Utr/Nk B T m 


S tr /Nk B 


Na 
Cu 


0.65 
1.00 


0.415 
0.332 


0.72 
0.80 



Finally for the liquid, we must evaluate the transit con- 
tributions to the internal energy and entropy at melt. It 
was just this requirement, in the present application of 
DFT to liquid dynamics theory, that motivated our de- 
velopment of an improved model for the transit free en- 
ergy of monatomic liquids. This model consists of two 
parts: (a) The available high-T experimental entropy 
data were analyzed in terms of the entropy formulas, 
Eqs. © and (|A10[) . revealing a scaled T-dependence of 
S l tl .(V, T) at fixed volume [lj], and (b) a statistical me- 
chanical model for F^{V,T) was calibrated to this ex- 
perimental S l tT (V,T) function, yielding model equations 
for all thermodynamic functions which derive from the 
transit free energy |l2| . The model provides universal 
curves for S[JNk B and U l t jNk B T in terms of T/6 tI {V), 
where 9ti{V) is the material-specific scaling temperature 
for the transit entropy. Results for Na and Cu at melt 
are listed in Table |TJ Volume dependence of the transit 
contribution is neglected. 

The temperature 9ti{V) plays a role in the transit con- 
tribution similar to that played by 0o(V) in the vibra- 
tional contribution; it sets a material-specific tempera- 
ture scale. We know how to calculate 0o(V) from first 
principles because the relevant term in the liquid Hamil- 
tonian (the vibrational part) is well-understood; however, 
the theory for the transit contribution to the Hamiltonian 
is still under development. Once the transit term is un- 
derstood, we hope to be able to calculate 9tr(V) from 
first principles as well. For now, parameterization from 
data will suffice. 

In the crystal, the system moves in the crystal potential 
energy valley, and this motion is well described by lattice 
dynamics theory |l|. In the harmonic approximation, 
and neglecting valley-valley intersections, the equations 
of Appendix [5] apply. Therefore, the total crystal free 
energy is 

F C (V, T) = %{V) + T v c ib (V, T) + F^(V, T), (4) 

where $q(V) i s the system potential at the valley min- 
imum, the crystal structure, and T° ib (V, T) is the con- 
tribution from lattice vibrations. Again the electronic 
excitation term F^(V,T) is very small and is given to 
sufficient accuracy by free electron theory. 

Corresponding to Eq. Q, the crystal internal energy 
and entropy are given by 

U C {V,T) = %{V) + U^ h {V,T) + U^{V,T), (5) 

S C (V,T) = S c vib (V,T) + S c cl (V,T). (6) 



TABLE II: Setup parameters for the VASP calculations. The 
Monkhorst-Pack k-mesh is recorded as [n, n, n], followed by 
the number of k-points in the irreducible Brillouin zone as 
(%). The quantity "translational invariance" is the maxi- 
mum magnitude of translational eigenvalues relative to that 
of the lowest pure vibrational mode. "DM" is short for "dy- 
namical matrix." 



Quantity 



Na 



Cu 



valence electrons 1 11 

planewave cutoff [eV] 101.7 341.6 

max core radius [A] 2.5 2.3 

EDIFF [eV] 10~ 8 10" 8 

k-mesh for E l (V) [3, 3, 3] (14) [5, 5, 5] (63) 

k-mesh for liquid DM [3, 3, 3] (14) [2, 2, 2] (4) 

translational invariance (liquid) 10 -6 10 -6 

crystal structure bcc fee 

k-mesh for crystal DM [3, 3, 3] (10) [2, 2, 2] (2) 

translational invariance (crystal) 10 -13 10 -7 



For Na and Cu, as with most elements, the high-T ex- 
pansions Eqs. (|A10[) and (|Allj) are accurate at T m for the 
crystal as well as the liquid. The fundamental differences 
in finite- N errors for crystal and liquid are described in 
the following Section. 



III. ELECTRONIC STRUCTURE 
CALCULATIONS 

A. Calculations for the liquid 

Our supercell consists of 150 atoms in a cubic box with 
periodic boundary conditions. The value N = 150 is 
large enough that finite-V errors are not serious, and 
small enough that a sufficient number of total energy 
calculations (a few thousand) can be done for each ele- 
ment. 

The DFT calculations are done with the VASP code 
, using the projector augmented wave (PAW) method 
in the generalized gradient approximation (GGA) 
The k-point mesh was automatically generated us- 
ing the method of Monkhorst and Pack [l8[ . The setup 
parameters are listed in Table [TTJ It is the large size of 
the real-space supercell which allows us to use few k- 
points in comparison to the large number (several thou- 
sands) needed for crystal metal calculations with small 
unit cells. 

To locate structures, we followed the procedure out- 
lined in Sec. HI [lfl, [Il(. A quench is considered con- 
verged when the energy decrease in one iteration is less 
than 10 times the current EDIFF setting, where EDIFF 
defines the energy convergence criterion used in the SCF 
procedure (see Table HI]) . Initially, quenches were done 
from a separate stochastic configuration at each volume. 
To save computer time we quenched to one structure, 



FIG. 1: For Na at N = 150: DFT results for the structural 
energy E versus volume V for the liquid (E l ), and for the 
crystal (E c ). 
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FIG. 2: For Cu at N = 150: DFT results for the structural 
energy E versus volume V for the liquid {E l ), and for the 
crystal (E c ). 
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then scaled this structure to slightly larger and smaller 
volumes, and quenched these configurations to new struc- 
tures. This was repeated until the desired range of vol- 
umes was covered. In quenching with DFT, the follow- 
ing procedure saves computer time: quench to conver- 
gence at Monkhorst-Pack k-mesh [1, 1, 1], then take 
Monkhorst-Pack k-mcsh [2, 2, 2] and quench to conver- 
gence again, and so on. This is faster, not because it re- 
quires fewer quench steps, but because most of the steps 
are at a smaller number of k-points. 

The DFT energy of a random structure is denoted 
E l (V). The results for Na and Cu are shown in Figs. Q] 
and [2] respectively. The crystal structure energies E C {V) 
are also shown in the figures. From the potential energy 
hypothesis mentioned in Sec. HH for a given element at 
a given V, the random structural energies should occupy 
a distribution whose width is small compared to fcsT m , 
and whose mean lies above E c by around ksTm. This 
characteristic is well verified for Na at the volume V^ of 
the liquid at melt US, EH S3] • Here, since each E l {V) is a 
representative of the random distribution at that V , the 
curves of E l (V) and E C (V) in Fig. Q] confirm the char- 
acteristic random structure distribution over a range of 
volumes for Na. The same confirmation is provided for 
Cu in Fig. [2] These energies are then normalized as de- 
scribed in Appendix iBlto provide the structural potential 
for the liquid, $|,(V), 

In the many quenches done in the present study, very 
few symmetric structures have appeared. These struc- 
tures invariably have DFT energies lying noticeably be- 
low the curve of E l (V), and above E C (V). This also ac- 
cords with the potential energy hypothesis of V-T theory, 
Sec. lU and accords with previous findings [Hi LL9|, [20] for 
Na at V^. 

For a given liquid at fixed density, the macroscopic ran- 
dom structure is better approximated with an ever larger 
supercell. Hence at finite N all calculated potential en- 
ergy parameters will have an error due to finite resolution 



of the structure. Moreover the vibrational characteris- 
tic temperatures are subject to a second finite- A^ error, 
due to the incomplete resolution of the frequency distri- 
bution. This second error is the dominant error in our 
liquid characteristic temperature calculations. 

Among the normal modes are three, A = 1, 2, 3, repre- 
senting uniform translation of the system. Their eigen- 
values are in principle zero due to translational invari- 
ance (also called the acoustic sum rule in crystal theory). 
In practice the translational eigenvalues are zero only to 
numerical accuracy, and are of either sign. The remain- 
ing 3A" — 3 modes are pure vibrational and have positive 
eigenvalues, by the definition of a structure as a local 
minimum. To check translational invariance, the ratio of 
eigenvalues \w%\ /ui\ is calculated for A = 1,2,3, where 
mode A = 4 is the lowest lying pure vibrational mode. 
The maximum value of this ratio in our final calculations 
for each element is listed under the designation "transla- 
tional invariance" in Table [TTJ The requirement is clearly 
satisfied to high numerical accuracy. 

Here we are interested in moments of the frequency 
distribution, which are related to characteristic temper- 
atures 9 n by Eqs. (|A3|) - (IA5|) . where the translational 
eigenvalues are excluded from the average. Our calcula- 
tions of these 9 n for Na are shown in Fig. [31 Also shown 
are values from a well-tested interatomic potential for 
Na, for a 500 atom system at the volume of the liquid at 
melt |2(|. Agreement is excellent for 9q, and is good for 
9\ and 9^. It is common for theory and experiment both 
to give more reliable values for 9$ than for other 9 n . The 
reason is that #o, being the log moment, is uniformly sen- 
sitive to all oj\ in the spectrum, while other 9 n are more 
sensitive to higher (or lower) frequencies, hence depend 
on a smaller portion of the spectrum. 



FIG. 3: For liquid Na at JV = 150: DFT results for the 
vibrational characteristic temperatures 6 n versus volume V, 
for n = 0,1,2 (open symbols). Also the same 6 n at TV = 
500 and V = V^, from Na interatomic potentials [2(| (solid 
symbols) . 



TABLE III: Intermediate Theoretical Results 
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B. Calculations for the crystal 



Quantity 



Crystal at T = 



V^ [A 3 /atom] 
$o [meV/atom] 
U vih [meV/atom] 
U\ T [meV/atom] 
U' cl [meV/atom] 
'S'vib [fcs/atom] 
Sl T [fes/atom] 
S' el [fcs/atom] 



Na 



40.93 

4.62 

96.65 

13.27 

0.82 

6.96 

0.72 

0.05 



Cu 



V T c Bf [A 3 /atom] 


37.24 


12.02 


D [eV/atom] 


1.28763 


3.7047 


|/Vfc s 6»i [meV/atom] 


16.05 


29.3 


ffc s #i/3fc B T m 


0.167 


0.087 


Crystal at T m 






V£ [A 3 /atom] 


39.79 


13.05 


$o [meV/atom] 


-12.22 


3.1 


Uvih [meV/atom] 


96.72 


351.7 


C/g! [meV/atom] 


0.82 


5.3 


•S'vib [fcs/atom] 


6.79 


9.03 


5J [fes/atom] 


0.05 


0.09 


Liquid at T m 







13.54 

96.5 

351.7 

38.9 

5.2 

9.39 

0.80 

0.09 



For the structural energy calculations, we use one 
primitive unit cell with periodic boundary conditions, 
and increase the k-mesh to convergence. This represents 
the infinite lattice model [lj, and there are no finite- TV er- 
rors. The situation contrasts with the liquid calculation, 
where the structure itself contains a finite- N error. 

For the crystal, our calculated phonon moments can 
be tested against results from inelastic neutron scatter- 
ing 21] . This will ultimately allow us to estimate the 
accuracy of the calculated liquid moments. To this end, 
we calculate the crystal dynamical matrix for a large su- 
percell with periodic boundary conditions, just as for the 
liquid. The procedure yields phonons with wavevectors 
commensurate with the supercell. The crystal structure 
is precise, but the phonon moments have finite- N error 
due to the limited resolution of the frequency distribu- 
tion. Since this is also the major error in the liquid, we 
expect the total error in vibrational characteristic tem- 
peratures to be about the same for crystal and liquid at 
a given N. 

To calculate the dynamical matrix, a rhombohedral 
supercell of 5 x 5 x 5 primitive unit cells was constructed 
for bcc and fee lattices [22[. This gives N — 125, as close 
as possible to the liquid value of N. A fixed k-mesh was 
chosen, the same as that used for the liquid; the same 
k-mesh produces a smaller number of k-points in the 
irreducible Brillouin zone for the crystal (see Table [TTJ) . 



IV. RESULTS 
A. Intermediate theoretical results 

For the crystal at T = and T m , and for the liquid at 
T m , we fit the four- parameter Vinet-Rose function [23[ 
to our calculated F(V,T) versus V. From this we find 
the volume at P — 0, and the isothermal bulk modulus 
B at that volume. Intermediate theoretical results calcu- 
lated from the formulas of Sec. [IT] are listed in Table [TTTT 
These show the relative importance of various theoret- 
ical contributions to internal energy and entropy. The 
discussion here is for Na and Cu collectively, and is qual- 
itatively applicable to monatomic crystals and liquids in 
general. 

At T = 0, the free energy is given by Eqs. (|Blj) - (|B3l) . 
For accurate theoretical work, the zero-point energy can- 
not be neglected. For the light elements the zero-point 
energy measurably affects the volume at P — (see also 
Ref. [24j, Table 16.3). For all elements, the zero-point 
energy is important in the internal energies of crystal 
and liquid states. This is shown by the ratio of the zero- 
point energy to the classical vibrational energy at melt, 
|fcs^i/3fcsT' m , listed in Table ILTT1 If the zero-point en- 
ergy is omitted from theory, this ratio is the relative er- 
ror made in the (dominant) internal energy contribution 
U vih {V,T). 

For the crystal at melt, the internal energy and entropy 
are given by Eqs. ([5]) and ©. From Table [TTTl the domi- 



nant energy contribution is U^ ih . The small contribution 
from $3 is a combination of the volume-dependent part of 
E C (V), and the zero-point energy, according to Eqs. (|B2j) 
and ([53]l . The contribution £? C (F) -£ c (V; c cf ) can be read 
from Figs. [T] and [2] The dominant entropy contribution 
is again vibrational. At T > T m , S^ ib depends almost 
entirely on T/0§(V), from Eq. (jMOj. In Table Hill the 
electronic contributions are quite small, being < 2% for 
the crystal at melt. These contributions are much larger, 
say up to 10%, for metals with unfilled d-bands 1251 ]. 

The liquid thermodynamic functions contain terms 
analogous to those in the crystal, plus an added con- 
tribution from transits, Eqs. ^Q - ^ and (|B4l) . For the 
liquid at melt, the character of contributions from the 
structural potential $q, from vibrations, and from elec- 
tronic excitations is qualitatively the same as described 
above for the crystal at melt. Again the contribution 
from E l (V) - E c {V r c c{ ) can be read from Figs. [Q and d 
However, the transit contribution in Table ITO1 is approx- 
imately 10% for the liquid at melt, and is therefore im- 
portant for an accurate theory. This is the only entry 
in Table IIIII not obtained from electronic structure cal- 
culations. But this term also will be amenable to DFT 
calculation, as soon as a model for the transit Hamilto- 
nian is developed. 



B. Comparison of theory and experiment 

Comparison of theory and experiment for the crystal 
at T = and P — is listed in Table IIVI Differences are 
expressed in the quantity A, defined in general by 



A 



theory — expt 
expt 



(7) 



For the crystal volume and bulk modulus, the agreement 
is excellent, at the customary level for ab initio crystal 
calculations (see e.g. Ref. [2J], Tables 16.2 and 16.3). 

In Table HVl comparison of d^(V) for n = 0, 1, 2 is at 
the experimental volume V£ eas (see Ref. [24]], Table 15.1). 
The experimental error in #£(V^ eas ) is estimated to be 
0.1 - 0.5% (p. 151 of Ref. [Hj]). In our experience, lattice 
dynamics theory in ab initio evaluation can account for 
the experimental B c n to an accuracy around 1% at best. 
To achieve such accuracy is not our goal here. Accuracy 
of the ab initio 9^ in Table HVl is quite respectable, with 
A in the range —2% to +5%. Error at this level is a 
minor effect in the comparison of theory and experiment 
for the crystal at melt. Moreover, we attribute the theo- 
retical error mainly to small system size (TV = 125), since 
the corresponding small number of vibrational modes can 
only poorly represent the actual crystal frequency distri- 
bution. This problem is easily remedied by increasing 
N. 

Comparison of theory and experiment for the crystal 
at melt, and for the liquid at melt, is listed in Table fVl 
Notice the nuclear motion causes thermal expansion, as 
seen in the volume at melt. Specifically, V^ is larger than 



V at the minimum of E C (V), and V^ is larger than V 
at the minimum of E L (V), Figs. [T]and[5] The volumes 
are in excellent agreement with experiment. Also, at this 
point we can see that the volume errors are systematic: 
for all three states, the crystal at T = and at T m , and 
the liquid at T m , the volume error is —0.01 for Na, and 
+0.03 for Cu. The bulk modulus, being essentially the 
curvature of $q(V), for crystal or for liquid, has larger 
error than the volume itself. Moreover, the experimental 
and theoretical determinations of B in Table [V] are both 
subject to significant errors. The agreement of theory 
and experiment for B is as good as we can expect at T m . 

It remains to discuss the comparison of theory and ex- 
periment for the internal energy and entropy at melt, 
Table [V] The errors here are sufficiently small that in- 
dividual A values cannot be interpreted. The mean and 
standard deviation of A for energy and entropy for crys- 
tal and liquid states is A = —0.005 ± 0.019. Hence the 
errors are essentially pure scatter. Contributions to the 
scatter arise from two major sources, experimental error 
at the level of 0.005 - 0.010, and computational error due 
to small system size at the level of 0.01 - 0.02. Addi- 
tional smaller errors result from the slightly inaccurate 
theoretical volume, from error in the free-electron model 
for electronic excitation, and from neglect of electronic 
excitation- nuclear motion interaction [3ll . l32j . It follows 
that the results in Table [V] for internal energy and en- 
tropy at melt are consistent with known errors. 

There is one more systematic property of the compari- 
son in Table IVl that holds separately for each metal. For 
the bulk modulus, A is large and positive and is roughly 
the same for crystal and liquid, while for each remaining 
property, A is small and approximately the same magni- 
tude for crystal and liquid. The implication is that the 
liquid theory we study has the same level of accuracy as 
does lattice dynamics theory for crystals. 



V. CONCLUSIONS 

A. The present application of DFT to liquid 
dynamics 

Our primary conclusion from the present study is: Us- 
ing a standard implementation of DFT, as provided by 
the VASP code, it is possible to make ab initio calcula- 
tions of thermodynamic properties of monatomic liquids. 
The liquid calculations achieve the same level of accu- 
racy as does ab initio lattice dynamics for crystals, the 
most accurate crystal theory available. This conclusion 
is validated for Na and Cu, based on the comparisons of 
theory with experiment listed in Table |V] The following 
discussion adds detail to the primary conclusion. 

The comparisons for the liquid in Table [V] are at zero 
pressure and temperature T m . However, we fully ex- 
pect the agreement of theory and experiment to hold for 
a wide range of pressures and temperatures, or equiva- 
lently, a wide range of volumes and temperatures. The 



TABLE IV: Comparison of theory and experiment for the crystal at low T. Vf et and B c are at T = and P = 0. 0„ are at 



1 in 



the volume of experimental measurement. 



Quantity 


Theory 


Na (bcc) 
Expt. 


A 




Theory 


Cu (fee) 
Expt. 




A 


V T c e{ [A 3 /atom] 


37.24 


37.68° 


-0.012 




12.02 




11.70" 




0.027 


B c [GPa] 


7.76 


7.3-7.6' 


0.042 




138.7 




142.0 C 




-0.023 


KSeas [A 3 /atom] 


37.98 


37.98" 


— 




11.70 




11.70" 




— 


<?o [K] 


111.2 


113.3 d 


-0.019 




236.5 




225.3 d 




0.050 


0i [K] 


161.7 


163 d 


-0.008 




330.2 




315 d 




0.048 


*S [K] 


165.8 


166 d 


-0.001 




332.5 




317 d 




0.049 


"Tables 15.1 and 19.1 of Ref. [H| 


















6 See Refs.Jimi 




















c See Ref. 
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d See Ref. 


21] and Tables 15.1 and 19.1 of Ref. J2J 


















TABLE V: 


Comparison of theory and experiment for 


crystal 


and liqu 


id at T„ 






Quantity 


Theory 


Na (bcc) 
Expt. 


A 


Theory 


Cu (fee) 
Expt. 




A 




Crystal Data 




















V c [A 3 /atom] 


39.79 


40.27" 


-0.012 


13.05 




12.62" 




0.034 




B c [GPa] 


6.47 


5.8 6 


0.116 


97 




90 c 




0.078 




U c [meV/atom] 


85.32 


89. l d 


-0.042 


353.1 




358.9 d 




-0.016 




S c [fcs/atom] 


6.84 


6.93 d 


-0.013 


9.12 




8.93 d 




0.021 




Liquid Data 




















V' [A 3 /atom] 


40.93 


41.27" 


-0.008 


13.54 




13.19" 




0.027 




B ! [GPa] 


5.93 


5.3 C 


0.119 


85.7 




73.6 / 




0.164 




[/' [meV/atom] 


115.36 


115. 9 d 


-0.005 


492.3 




492. 7 d 




-0.001 




S 1 ' [fcs/atom] 


7.73 


7.78 d 


-0.006 


10.28 




10.09 d 




0.019 





"Tables 19.1 and 21.1 of Ref. 

6 See Ref. [11 

c See Ref. M 

d See Ref. H 

e See Refs.jpllfl 

•''See Ref. SI 



volume dependence is contained in functions calculated 
by DFT, primarily & (V) and l Q (V). We can expect 
these functions to be as accurate for the compressed liq- 
uid as for the liquid at zero pressure. Moreover, at the 
fixed volume Vj^ , the temperature dependence of the ex- 
perimental thermodynamic data is accurately accounted 
for by the equations of Sec. |TT] and Appendix [A] This 
is because those equations, and the experimental data 
for entropy at high temperatures, were used to calibrate 
the statistical mechanical model for the transit free en- 



ergy [14J. Hence the level of agreement between theory 
and experiment found in Table IVl should persist to higher 
pressures and temperatures. 

A potentially useful comparison can be made between 
the present technique and ab initio MD. At a given N, ab 
initio MD requires far greater computer resources to cal- 
culate a similar set of thermodynamic data. Or, with a 
fixed computer resource, the present technique can study 



a larger system, and hence obtain greater accuracy by re- 
ducing finite- TV errors. On the other hand, ab initio MD 
data contain both vibrational and transit contributions. 
Both techniques together can separate the transit and 
vibrational contributions to statistical mechanical func- 
tions. The combined techniques have the potential to 
reveal the physical nature of transit motion. 

What is achieved by making DFT calculations for the 
crystal at T = 0? First, before we can compare theory 
and experiment for any condensed matter state, it is nec- 
essary to adjust the DFT energy calculations to have the 
thermodynamic zero of energy. This is accomplished by 
means of Eqs. (|B2I) and (IB3|) . A second useful result is 
the confirmation that DFT calculations are accurate for 
$o(F), and for the characteristic temperatures 6q, Q\, 



and #2 • This is shown by the comparisons of Table IIVI 
This confirmation lends confidence to similar calculations 
for the liquid. 
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What is achieved by comparing theory and experiment 
for the crystal at T m l This comparison is rarely per- 
formed in crystal physics research, and is of interest in 
itself. The comparisons of Table IVl confirm that the DFT 
calculations are accurate for ^(V) and 6q(V), and their 
volume derivatives. This confirmation shows the level of 
accuracy which ab initio lattice dynamics can achieve for 
the crystal at melt, and it also lends confidence to the 
similar calculations for the liquid. 



B. Verification of V-T theory 

The primary conclusion is based on the comparison of 
theoretical and experimental data, and uses no condi- 
tion on the validity of theory. But the well-established 
reliability of DFT calculations supports a secondary con- 
clusion regarding the theory itself. The theory consists 
of two parts, the V-T equations and DFT calculations. 
We can reasonably assume that DFT is as accurate for 
the liquid as for the crystal. Then, at this level of accu- 
racy, the comparisons of theory and experiment for the 
liquid, in Table |V1 confirm the equations of V-T theory 
as described in Sec. HU Moreover, confirmation of the 
equations provides confirmation of the nuclear motion de- 
scribed by the equations. The following discussion adds 
detail to this secondary conclusion. 

In the formula for internal energy, Eq. ©, the major 
contribution is $ l (V) + U l vih {V, T). This is a purely the- 
oretical function, and describes normal-mode vibrational 
motion of the nuclei. Also in Eq. ([5]), no significant error 
is contributed by the small terms £//,.(¥, T) and U^(V, T). 
Agreement of theory and experiment for the internal en- 
ergy for one liquid demonstrates consistency of the vi- 
brational motion with experiment. The same argument 
applies to the liquid entropy, Eq. © , where the same nu- 
clear motion is described by the purely theoretical func- 
tion S vih (V, T). Agreement of theory and experiment for 
the entropy for one liquid demonstrates consistency of the 
vibrational motion with experiment. Further, the quan- 
tities calculated by DFT, $p(V) for the internal energy 
and l (V) for the entropy, are not theoretically related, 
so the confirmations obtained from internal energy and 
entropy are independent. 

This situation is analogous when one compares theory 
with experiment for the liquid volume and bulk modulus. 
The theoretical functions tested are volume derivatives of 
$o(V) and (V), and the comparison with experiment 
for Vj^ and i3(V^,T m ) are independent tests. 

Altogether in Table [V] four independent consistency 
tests of the nuclear motion are provided for each of two 
elemental liquids. The agreement of theory and experi- 
ment within small errors for all these tests provides the 
following two-part verification of V-T theory (for Na and 
Cu): (a) Many-body harmonic vibrational motion de- 
scribed by the parameters $ (V) and 8 l (V) is the dom- 
inant contribution to the thermodynamic functions, and 
(b) While the experimental liquid moves rapidly among 



all the valleys in the potential energy surface, a single 
random valley at each volume serves to calculate the vi- 
brational motion and its contribution to thermodynam- 
ics. 



Appendix A: Statistical mechanics for a single 
potential energy valley 

The system has N atoms in a volume V at temperature 
T, and is confined to a single harmonic potential valley. 
The system potential energy at the valley minimum, the 
structure, is <&o(V). The normal mode vibrational fre- 
quencies are uj\(V) for A = 4, ... , 3JV, where the three 
translational modes are omitted from statistical mechan- 
ics. The Hclmholtz free energy is F, given by 



F(V,T) 

F vih (V,T) 



E 



\hu x 



F vih (V,T), 
k B T\n(n\ 



(AI) 
,(A2) 



where n\ is the Bose-Einstein distribution function. In 
low- and high-T regimes, i 7 ^ depends on only a few 
characteristic temperatures 6 n , which are related to mo- 
ments of the frequency distribution. Here the important 
9 n are given by 



\n(k B 8 Q ) 
k B 0i 



MM) , 

4 

o (M > 



(A3) 
(A4) 



k B 2 



(M ;, (A5) 

) is over the vibrational frequen- 



where the average 
cies. 

Functions useful for comparing with experimental data 
are the internal energy U(V, T), and the entropy S(V, T). 
These are obtained from the free energy, and Eq. (|A1[) 
yields 

S(V,T) = S vih (V,T), (A6) 

U(V,T) = $ (V) + [/ vib (V,T). (A7) 

We shall be interested in the vibrational contributions 
only at T = and T m . At T = 0, 



S vih (V,T = Q) = 0, 
U vih {V,T = Q) = lNk B 8i(Y) 



(A8) 
(A9) 



The right side of Eq. (JA9I) is just the harmonic zero-point 
energy. This equation will be needed below to evaluate 
the thermodynamic zero of energy. For most monatomic 
crystals and liquids, the nuclear motion is nearly classical 
at T > T m , and a high-T expansion of n\ in Eq. (IA2|) is 
valid. This expansion gives 



S vib {V,T) = Wk B 



In 



T 



h(V) 



1 



1 

40 



U vih (V,T) = 3Nk B T 



200 

T 



1 
20 



e 2 (v) 

T 



(A10) 



(All) 



If V c f is the crystal volume at T = and P = 0, then by 



The series starting with (QijT^) expresses the quantum 
corrections, and only this leading term is required in the 
present study. 



Appendix B: Normalizing DFT energies for 
thermodynamics 

To evaluate the functions in Appendix [Al we require 
the potential energy at the valley minimum <&q(V) and 
the moments of the vibrational spectrum 8 n for n — 
0, 1, 2. Extraction of the moments from DFT is described 
in Sec. HU <&$(V) is determined as follows. 

All system energies are to be measured from the ther- 
modynamic zero of energy, which is the energy of the 
crystal at T — and P = 0. Here P is pressure, given by 
P = — (dF/dV) T . At T = 0, the crystal free energy is 



f c (v) = u c {v) = sg(V) + -Nk B ei(y). 



(Bl) 



DFT calculations provide total energies E(V) measured 
with respect to isolated atoms. To correct the E(V) to 
the thermodynamic energy zero, we denote the DFT crys- 
tal structure energy by E C (V), and define the constant 
Dby 



definition U c (V r c cl ,T = 0) = 0. Equations (JBTJ) and (B2 
can be solved for D to find 



D 



£ c (K C ef) + z Nk B0i(V r c e{ ) 



(B3) 



Equations (|B2[) and (JB3I) define $o(^) m terms of quan- 
tities calculated by DFT. To observe the thermodynamic 
energy zero for a given atomic system, the same D is 
added to every energy calculated in DFT, for every con- 
densed matter phase, at all V. With this, only DFT 
energy differences enter the quantities we need to calcu- 
late. 

As for the crystal valley, potential energy properties 
of random valleys depend on the system volume. Hence 
to include volume dependence, a representative random 
valley is required at a range of volumes. Equation (|B2[) 
holds for the liquid, i.e. for each random valley in the 
form 



& (V) =E\V) + D. 
Acknowledgments 



(B4) 



$g(V) =E C (V) + D. 



(B2) 
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